Exact Diagonalization Approach for the D = oo Hubbard Model 
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We present a powerful method for calculating the thermodynamic prop- 
erties of the Hubbard model in infinite dimensions, using an exact diagonal- 
ization of an Anderson model with a finite number of sites. At finite temper- 
atures, the explicit diagonalization of the Anderson Hamiltonian allows the 
calculation of Green's functions with a resolution far superior to that of Quan- 
tum Monte Carlo calculations. At zero temperature, the Lanczos method is 
used and yields the essentially exact zero-temperature solution of the model, 
except in a region of very small frequencies. Numerical results for the half- 
filled case in the paramagnetic phase (quasi-particle weight, self-energy, and 
also real-frequency spectral densities) are presented. 

PACS numbers: 71.10+x,75.10 Lp, 71.45 Lr, 75.30 Fv 
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Following the pioneering work of Metzner and Vollhardt JT| , the limit of large dimensions 
for models of strongly correlated fermions has received much attention. In this limit, the 
highly intricate quantum many-body problem simplifies considerably and leads to a non- 
trivial mean- field theory [0. Remarkably, this limit captures many features of the physics 
in finite dimensions and gives a very successful description of quantum fluctuations. 

In spite of the considerable simplification obtained in taking the large D limit, the mean- 
field equations still have to be solved numerically. Up to now, all calculations 0,0,0 have 
relied on the Hirsch-Fye Quantum Monte Carlo (QMC) algorithm |J. A major limitation of 
this scheme is the difficulty of accessing the low-temperature regime, where statistical and 
finite time-step discretization errors of the QMC algorithm become very important. 

In this paper, we present a powerful method for solving these mean-field equations, which 
leads to an essentially exact solution in the imaginary frequency domain. As an example, we 
consider the Hubbard model on a lattice of infinite connectivity z — > oo which, after proper 
rescaling of the kinetic energy, is written as 

H = - ~7^= c t c o + h - c - + U Y, n ^ n iU (!) 

<ij>a V ^ Z i 

The calculation of the single-site properties of the Hubbard model in this limit reduces to 
the self-consistent determination of the on-site Green's function G(u>) of the Hubbard model 
and of a bath Green's function Gq(u), which describes the interaction on the single site 
with the external environment. G(uS) and Gq(u) are related by a self-consistency condition 
which, on the Bethe lattice, reads: 

G Q \u) = uj + fi- 1 -G(u) (2) 

It is for simplicity only that we restrict our attention in this paper to the z — > oo Bethe 
lattice. 

As is well known R] ||, the on-site Green's function of the Hubbard model may be 
interpreted as the Green's function of an Anderson model 

H And = e d J2 d a d <r + e * a tor a *w + Un d^ dl + ( V k a ka d <r + h -C.) (3) 

a a,k=2 a,k=2 

in which the function Go(iuj n ) is given by the U — Green's function of the impurity 

n s y2 

G Q (iuj n ) = G^"Vn) = K - e d - ii - £ — ^ (4) 

Given the infinite number of degrees of freedom of the models defined in eq. (|l|) and eq. ([|), 
it is evident that strict self-consistency can only be obtained with a continuous Anderson 
model, i. e. with n s = oo. The main result of the present paper is that a systematic 
approximation of Go(iu) with a finite-n s Anderson model gives extremely good results. We 
stress from the beginning that we are interested in an approximation of the imaginary- 
frequency Green's functions only. 

In practice, we approximate any Gq 1 (iuj) by a function Gq 1 And (iu>) with a finite number 
n s of sites. This can be cast into a minimization problem in the variables and For 
this paper, we choose the following cost function: 
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' - J2 \G Q l {i^ n ) - G Q lAnd {iuj n )\ z (5) 

-L n 



Umax x n =0 



where n max is chosen sufficiently large {oJ nmax » maxk(ek)) |§. We search for the param- 
eters efc and Vk minimizing the \ 2 in eq. (|5|) with a standard conjugate gradient method 

For a small number of sites, n s < 6, the Green's function G(iu n ) can be obtained exactly 
from the complete set of eigenvectors and eigenvalues of the Anderson Hamiltonian eq. @. 
The procedure 

G-^ n ) ^ G 1 A ^{iu n ) eq -M G{iu n ) eq -M G- Q \iu n ) (6) 

is then iterated to convergence. 

The following observations are made: 

1) We notice in general very small differences between G~ l (iuj) and G~ 1 d (iuf) as ex- 
pressed by small minimal values of x 2 m eq. ©• X 2 decreases by approximately a constant 
factor each time we add one more site. 

2) The extensive comparisons with QMC results || which we have undertaken indicate 
that exact diagonalization is by far the superior method for this problem. As an example, 
we show in fig. [I] QMC and exact diagonalization data for the half-filled Hubbard model at 
(3 = 32 and U = 3. The Monte Carlo data are shown for a imaginary-time discretization of 
At = 1, .5, and .25 (c/j e. g. |J), and the exact diagonalization data for n s = 3,5. It may 
be worthy of notice that the diagonalization calculations can be obtained in a few minutes 
on a work station, while for the QMC data acquisition (at At = 0.5) several hours were 
needed (several days for At = 0.25). 

Beyond n s = 6, the size of the Hilbert space becomes too large for an explicit diagonal- 
ization of the Anderson Hamiltonian. However, the calculation of zero-temperature Green's 



functions is still possible by means of the Lanczos algorithm ||11|| , which allows us to easily 
calculate G(ito) and G${iuj) up to n s ~ 10 on a workstation ||12|| . The fit with the An- 
derson model is performed as before. We simply replace the Matsubara frequencies by a 
fine grid of imaginary frequencies, which correspond to a "fictitious" inverse temperature 
(3 {oj n = (2n + 1)tt/P). (3 introduces a low-frequency cutoff in an obvious way. In fig. || 
we display the functions G^iiu) and G~ 1 And {iu) for U = 2 and U = 4.8. At the scale 
of the figure the two curves can hardly be distinguished, and an essentially perfect fit (i.e. 
perfect self-consistency) is obtained in the whole range of frequencies. The inset in the figure 
shows a blow-up of the small frequency regime as the number of sites is increased. Notice 
the systematic amelioration of the fit. Furthermore, the 'physical' Green's function Gq 1 is 
extremely independent of n s , especially at high frequency. Already at uj = 0.11, e.g., Gq 1 
varies by less than 0.0001 between n s = 6, 8, and 10. 

For the data at U = 4.8, the quality of the fit is excellent even with a small number 
of sites. This is easily explained by the existence of a physical cutoff in frequency, which 
results from the Mott gap. 

We now pass to the calculation of other physical quantities and present in fig. ^ some 
results for the quasi-particle spectral weight Z calculated from the slope of the self-energy 
S = Gq 1 — G^ 1 . In the inset of fig. |3] we present the raw data of ImT,(iu) at small frequencies 
from which the spectral weight is extracted (ImE(iuj) ~ (1 — 1/Z)u + ...). To get a truly 
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stabilized slope of £ we have found it to be necessary to reach very low temperatures. The 
main plot compares the results at n s = 10 with the "iterated perturbation theory" (IPT) 
result. This method is based essentially on the use of a weak coupling calculation to second- 
order in U of £ and has shown to give a satisfactory interpolation between the small and 



large U limits (exclusively at half filling and in the paramagnetic phase) |[L3| , ||14|| . On a 
few points we give in addition the results of the exact diagonalization at n s = 6 and n s = 8. 
Given the extremely good agreement between the values of Z calculated with n s = 8 and 
10, we are very confident of the numerical values presented. 



As discussed in ref [131, the IPT approximation leads to a first-order Mott-Hubbard 
transition (c/fig. |3|), and the quasi-particle weight Z jumps discontinuously at U ~ 3.6. We 
have only found limited evidence for such a scenario within the present approach. At n s = 6, 
we are unable to stabilize two solutions at the same values of the physical parameters (the 
coexistence of two solutions is indicative of a first-order phase transition). At n s = 8, and 
using a fictitious temperature of j3 = 120, we find a coexistence region within a very small 



interval of U: 4.45 < U < 4.60 ||15|| . Even though the question of the order of the transition 
will have to await a more detailed investigation, it seems to us to be difficult to reconcile 
our numerical results with a abrupt transition anywhere close to U — 3.6. 

Finally, we show some data concerning the one-particle spectral densities p(u>) = 
—ImG(uj + ie)/7c as obtained from the Lanczos calculation together with IPT-approximation 
solutions [|L|. Fig. [| shows the spectral density (n s =10) for different values of U. In the 
Fermi-liquid regime the spectrum of our finite-size Anderson model consists of a large num- 
ber of peaks, while in the insulating phase we systematically observe a simpler structure 
made of only a few peaks. As U is increased we see that p(uj) develops three well-separated 
structures: a central quasi-particle feature and two broad high-energy satellite features cor- 
responding to the formation of the upper Hubbard band. At U = 4.8 a gap is observed in 
good agreement with the approximate IPT solution. In the insets of Fig. ^ we also present 
the integrated single particle density of states corresponding to Lanczos and IPT solutions. 
The agreement between both curves is seen to be very good, provided we average over a 
frequency interval of uj ~ 0.5. This indicates that the calculated spectral density contains 
coarse-grained information about the exact solution, as it should be. Due to the discrete 
nature of the Anderson model used, the fine details of the spectrum are poorly reproduced. 

A remark is in order here: As is well known, the continuation of numerical data from the 
imaginary axis onto the real-frequency domain is a very difficult problem and constitutes 
for example one of the major limitations of the QMC method. Here we encountered the 
analytic continuation problem in the 'easy' direction. Indeed, in the present work very 
precise imaginary frequency data (c/fig, [TJ, fig. |^) can be obtained with a representation 
in uj , which very cleary has its limits (c/fig. [|). It is crucial in this context that we only 
attempt to satisfy the self-consistency condition on the imaginary axis. 

In conclusion, we have presented a powerful numerical method for simulating the D = oo 
models for strongly correlated fermions based on a self-consistent single-impurity model 
treated by exact diagonalization. At the temperatures reachable by quantum Monte Carlo 
calculations we get essentially the exact solution of the model. At lower temperatures 
unreachable by QMC, we get also an amazingly good solution, except in the region of very 
small frequencies where some difficulties appear due to the finite degrees of freedom in the 
representation of the free propagator (Gq 1 And ) of the impurity. 
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Elsewhere [|16[ we present a study of the instability of the normal phase with respect 
to superconductivity of an infinite-D two-band Hubbard model fTfll . That work and addi- 
tional calculations on the Hubbard model away from half-filling clearly show that the exact 
diagonalization method presented in this work is in no way limited to the particle-hole sym- 
metric point of the Hubbard model. Broken-symmetry phases, magnetic fields, etc, as well 
as the calculation of susceptibilities |TI| can be easily handled within this approach, which 
we expect to rapidly become a standard tool for the investigation of D = oo systems. 
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Figure Captions 



1. Comparison of G(r) for U — 3, and /3 = 32 between exact diagonalization (n s = 3,5 
(bottom solid lines) the results cannot be distinguished on the scale plotted) and 
Quantum Monte Carlo (solid: At = 1, dashed: At = 0.5). Inset: G(t = 4) vs. At 2 
(the QMC algorithm converges roughly in At 2 ). At At = 0. exact diagonalization 
results for n s = 3, 5. 

2. Plot of Gq 1 {ioj) and of G~ 1 And for two values of the interaction, U = 2 (Fermi liquid- 
regime) and U = 4.8 (Mott insulator regime) at n s = 10. The inset gives Go(^) -1 
and G~ 1 And at small frequency for n s = 6, 8, 10. Note the systematic improvement of 
the solution. The misfit between the two functions is the only source of error of the 
exact diagonalization approach. 

3. Quasi-particle weight Z as a function of U for the half-filled Hubbard model. The curve 
gives the IPT approximation, which predicts a first-order transition. The crosses give 
the results for n s = 10, with the corresponding results for n s = 6, 8 at two points. 
The inset shows the small-u; behavior of ImE(iuj) for n s = 6,8,10 from which the 
quasi-particle weight is calculated. Note the excellent convergence with n s . 

4. Density of states p(u) for different values of U (a value of e = 0.01 is used). We 



compare with the IPT density of states [13]. Insets: Comparison of the integrated 



densities of states between exact diagonalization and IPT. 
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